Introduction au Reservoir Computing avec ReservoirPy¶

Nathan Trouvain
Inria - Mnemosyne




PFIA - 28 juin 2021

 Sommaire¶

  • Concepts et fonctionnalités clés
  • Chapitre 1 : application à des données simples
  • Chapitre 2 : déployer les capacités génératives
  • Chapitre 3 : l'apprentissage en ligne
  • Chapitre 4 : démonstration d'application à des données réelles
  • Aller plus loin : comprendre les hyperparamètres

Concepts et fonctionnalités clés ¶

  • Numpy, Scipy, et seulement Numpy et Scipy
  • Exécution efficace (parallélisation des calculs, optimisation des règles d'apprentissage)
  • Plusieurs méthodes d'apprentissage, online, offline
  • Outils d'aide à l'optimisation des paramètres
  • Documentation: https://reservoirpy.readthedocs.io/en/latest/
  • GitHub: https://github.com/reservoirpy/reservoirpy

Informations générales¶

  • Le temps est toujours représenté par la première dimension des tableaux/vecteurs
  • Toutes les données (en entrée et en sortie) sont du type list[sequence], ou une séquence représente une série temporelle.

 Chapitre 1 : Application du Reservoir Computing a des données simples ¶

Série temporelle de Mackey-Glass

L'équation de Mackey-Glass décrit l'évolution de différent phénomènes biologiques, comme la quantité relative de globules rouges dans le temps.

Cette équation est une équation différentielle retardée :

$$ \frac{dP(t)}{dt} = \frac{a P(t - \tau)}{1 + P(t - \tau)^n} - bP(t) $$

où $a = 0.2$, $b = 0.1$, $n = 10$. Le retard $\tau$ est de $17$ pas de temps. $\tau$ a une forte influence sur le caractère chaotique de la série temporelle modelée. $17$ donne des résultats plutôt chaotiques.

In [3]:
from reservoirpy.datasets import mackey_glass

timesteps = 2510
tau = 17
X = mackey_glass(timesteps, tau=tau)

# rescale between -1 and 1
X = 2 * (X - X.min()) / (X.max() - X.min()) - 1
In [5]:
plot_mackey_glass(X, 500, tau)
  • Pas totalement imprédictible... (non aléatoire)
  • ...mais pas facilement prédictible non-plus (non périodique)
  • Comme la variation du rythme cardiaque, la météo, la bourse...

1.1. Tâche 1: prédire la valeur de la série à un horizon de 10 pas de temps¶

Prédire $P(t + 10)$ connaissant $P(t)$.

Préparation les données¶

In [7]:
from reservoirpy.datasets import to_forecasting

x, y = to_forecasting(X, forecast=10)
X_train1, y_train1 = x[:2000], y[:2000]
X_test1, y_test1 = x[2000:], y[2000:]

plot_train_test(X_train1, y_train1, X_test1, y_test1)

Préparation de l'ESN¶

In [8]:
units = 100
leak_rate = 0.3
spectral_radius = 1.25
input_scaling = 1.0
connectivity = 0.1
input_connectivity = 0.2
regularization = 1e-8
seed = 1234

In [10]:
Win = mat_gen.generate_input_weights(units, 1, input_scaling=input_scaling, 
                                     proba=input_connectivity, input_bias=True,
                                     seed=seed)

W = mat_gen.generate_internal_weights(units, sr=spectral_radius,
                              proba=connectivity, seed=seed)

reservoir = ESN(leak_rate, W, Win, ridge=regularization)

Entraînement de l'ESN¶

L'apprentissage est offline ("hors-ligne") : il n'a lieu qu'une seule fois, sur l'ensemble des données d'entraînement.

In [11]:
reservoir.train([X_train1], [y_train1], verbose=True)

Training on 1 inputs (2000 steps) -- wash: 0 steps

Train: 100%|██████████| 2000/2000 [00:00<00:00, 2537.93it/s]
In [12]:
reservoir.Wout is not None
Out[12]:
True

Test de l'ESN¶

In [13]:
y_pred1, states1 = reservoir.run([X_test1], verbose=True)

y_pred1 = y_pred1[0]

Running on 1 inputs (500 steps)

Run: 100%|██████████| 500/500 [00:00<00:00, 1439.26it/s]
In [14]:
plot_results(y_pred1, y_test1)

Coefficient de détermination $R^2$ et erreur quadratique normalisée :

In [15]:
r2_score(y_test1, y_pred1), nrmse(y_test1, y_pred1)
Out[15]:
(0.444170605212862, 0.23219740216886683)

 Amélioration des prédictions¶

  • La série de test est consécutive à la série d'entraînement.
  • Initialiser l'état du reservoir avec le dernier état d'entraînement devrait améliorer les performances.
In [16]:
states = reservoir.train([X_train1], [y_train1], return_states=True)

last_state = states[0][-1]
y_pred1, _ = reservoir.run([X_test1], init_state=last_state)

y_pred1 = y_pred1[0]
In [17]:
plot_results(y_pred1, y_test1)

Coefficient de détermination $R^2$ et erreur quadratique normalisée :

In [18]:
r2_score(y_test1, y_pred1), nrmse(y_test1, y_pred1)
Out[18]:
(0.9999562183846586, 0.003153349937792222)

 1.2 Compliquons la tâche¶

Passons d'un horizon de prédiction de 10 pas de temps à un horizon de 100 pas de temps

In [19]:
x, y = to_forecasting(X, forecast=100)
X_train2, y_train2 = x[:2000], y[:2000]
X_test2, y_test2 = x[2000:], y[2000:]

plot_train_test(X_train2, y_train2, X_test2, y_test2)
In [21]:
plot_results(y_pred2, y_test2, sample=400)

Determination coefficient $R^2$ and NRMSE:

In [22]:
r2_score(y_test2, y_pred2), nrmse(y_test2, y_pred2)
Out[22]:
(0.9686082619494742, 0.06591305950469496)

 Chapitre 2 : Déployer les capacités génératives des ESNs ¶

  • Entraînement de l'ESN sur une tâche de prédiction de 1 pas de temps.
  • Test de l'ESN sur ses propres prédictions (mode génératif).
In [23]:
units = 500
leak_rate = 0.3         # - leaking rate
spectral_radius = 0.99  # - rayon spectral
input_scaling = 1.0     # - facteur de mise à l'échelle des entrées (input scaling)
connectivity = 0.1      # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2  # et des entrées vers le reservoir
regularization = 1e-4   # - coefficient de régularisation (L2)
seed = 1234             # reproductibilité

 Entraînement à la prévision sur un horizon court¶

In [25]:
reservoir = reset_esn()

x, y = to_forecasting(X, forecast=1)
X_train3, y_train3 = x[:2000], y[:2000]
X_test3, y_test3 = x[2000:], y[2000:]

reservoir.train([X_train3], [y_train3])

 Génération¶

  • 100 pas de temps de la série de test utilisés comme "échauffement";
  • 300 pas de temps générés "à partir de rien";
In [26]:
start = 0; seed_timesteps = 100; nb_generations = 400

warming_inputs = X_test3[start:start+seed_timesteps]
X_t = X_test3[start+seed_timesteps: start+nb_generations+seed_timesteps]

X_gen, states, warming_out, warming_states = reservoir.generate(nb_generations, warming_inputs)
In [27]:
plot_generation(X_gen, X_t, nb_generations, warming_out=warming_out, 
                warming_inputs=warming_inputs, seed_timesteps=seed_timesteps)
In [28]:
states = reservoir.train([X_train3], [y_train3], return_states=True)

X_gen, states, _, _ = reservoir.generate(nb_generations, init_state=states[0][-1])

plot_generation(X_gen, X_test3[:nb_generations], nb_generations)

Chapitre 3 : Apprentissage hors-ligne ¶

Apprentissage se déroulant de manière incrémentale.

Utilisation de l'algorithme FORCE (Sussillo and Abott, 2009)

In [30]:
from reservoirpy import ESNOnline

Win = mat_gen.generate_input_weights(units, 1, input_scaling=input_scaling, 
                                     proba=input_connectivity, input_bias=True,
                                     seed=seed)

W = mat_gen.generate_internal_weights(units, sr=spectral_radius,
                              proba=connectivity, seed=seed)


reservoir_online = ESNOnline(leak_rate, W, Win, dim_out=1)  

reservoir_online.Wout is not None  # Wout est créée dès le début !
Out[30]:
True

Le dernier état calculé est cette fois-ci stocké dans l'objet. (Parallélisation impossible)

In [31]:
reservoir_online.state is not None
Out[31]:
True

Entraînement pas à pas¶

In [32]:
outputs_pre = np.zeros(X_train1.shape)
for t, (x, y) in enumerate(zip(X_train1, y_train1)): # pour chaque pas de temps de la série :
    
    # sortie avant entraînement, mise à jour de reservoir.state
    output_pre, state = reservoir_online.compute_output(x.reshape(1, -1)) 
    
    # entraînement à partir de reservoir.state
    reservoir_online.train_from_current_state(targeted_output=y.reshape(1, -1)) 
    
    
    outputs_pre[t, :] = output_pre
In [33]:
plot_results(outputs_pre, y_train1, sample=100)
In [34]:
plot_results(outputs_pre, y_train1, sample=500)

Entraînement sur une séquence complète¶

In [35]:
reservoir_online.reset_reservoir()

states = reservoir_online.train([X_train1], [y_train1])

pred_online, _ = reservoir_online.run([X_test1])  # Wout est maintenant figée
pred_online = pred_online[0]
In [36]:
plot_results(pred_online, y_test1, sample=500)

Determination coefficient $R^2$ and NRMSE:

In [37]:
r2_score(y_test1, pred_online), nrmse(y_test1, pred_online)
Out[37]:
(0.9999602106090204, 0.006511751721725012)

 Chapitre 4 : Application réelle : le chant des canaris domestiques ¶

Les données peuvent être téléchargées sur Zenodo : https://zenodo.org/record/4736597

In [38]:
im = plt.imread("./static/canary.png")
plt.figure(figsize=(5, 5)); plt.imshow(im); plt.axis('off'); plt.show()
In [40]:
display(audio)
Your browser does not support the audio element.

Plusieurs motifs temporels répétitifs differents à décoder : les phrases.

  • Un label par type de phrase à classifier dans le temps.
  • Un label SIL annotant le silence, à détecter également pour permettre la segmentation du chant.
In [41]:
im = plt.imread("./static/canary_outputs.png")
plt.figure(figsize=(15, 15)); plt.imshow(im); plt.axis('off'); plt.show()
In [44]:
X_train, y_train = X[:-10], Y[:-10]
X_test, y_test = X[-10:], Y[-10:]

 Entraînement de l'ESN¶

In [47]:
reservoir.train(X_train, y_train, verbose=True, workers=-1)

Training on 90 inputs (266054 steps) -- wash: 0 steps

Train: 100%|██████████| 266054/266054 [00:23<00:00, 11121.76it/s]
In [48]:
outputs, _ = reservoir.run(X_test, verbose=True)

Running on 10 inputs (32552 steps)

Run: 100%|██████████| 32552/32552 [00:02<00:00, 11873.05it/s]
In [50]:
scores  # pour chaque chant testé
Out[50]:
[0.9370629370629371,
 0.9093011305241521,
 0.9499136442141624,
 0.9412306983175847,
 0.962171052631579,
 0.9540372670807453,
 0.8858530661809351,
 0.9371994342291372,
 0.9174974217944311,
 0.8932743921692453]
In [51]:
print("Précision moyenne :", f"{np.mean(scores):.4f}", "±", f"{np.std(scores):.5f}")

Précision moyenne : 0.9288 ± 0.02471

 Pour aller plus loin¶

  • Essayer de varier le nombre de chants entraînés : combien faut-il de chants au minimum pour reproduire les résultats obtenus sur 100 chants ?
  • Trouver le meilleur jeu d'hyperparamètres : outils basés sur la bibliothèque hyperopt (https://github.com/hyperopt/hyperopt).
  • Essayer d'appliquer un input scaling différent pour chaque groupe de variables d'intérêt.

Merci de votre attention.¶

Nathan Trouvain
Inria - Mnemosyne




PFIA - 28 juin 2021

Aller plus loin : Comprendre les hyperparamètres et leurs effets ¶

In [52]:
units = 100             # - nombre de neurones dans le reservoir
leak_rate = 0.3         # - leaking rate
spectral_radius = 1.25  # - rayon spectral
input_scaling = 1.0     # - facteur de mise à l'échelle des entrées
connectivity = 0.1      # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2  # et des entrées vers le reservoir
regularization = 1e-8   # - coefficient de régularisation (L2)
seed = 1234             # reproductibilité

2.1. Le rayon spectral¶

Le rayon spectral est la valeur propre maximale de la matrice des poids du réservoir ($W$).

In [53]:
reservoir = reset_esn()

states = []
radii = [0.1, 1.25, 10.0]
for sr in radii:
    Win = mat_gen.generate_input_weights(units, 1, input_scaling=0.1, input_bias=True)
    W = mat_gen.generate_internal_weights(units, sr=sr, 
                                          proba=connectivity, seed=seed)
    reservoir.W = W
    reservoir.Win = Win
    s = reservoir.compute_all_states([X_test1[:500]])
    states.append(s[0])
In [54]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(radii)*100+10+i+1)
    plt.plot(s[:, :units_nb], alpha=0.6)
    plt.ylabel(f"$sr={radii[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()

  • $-$ rayon spectral $\rightarrow$ dynamiques stables

  • $+$ rayon spectral $\rightarrow$ dynamiques chaotiques

Rayon spectral et Echo State Property : rayon spectral $\rightarrow$ 1 (assure que les états internes ne sont pas affectés par l'initialisation).

2.2. Le facteur de mise à l'échelle des entrées (input scaling)¶

Il s'agit d'un coefficient appliqué à $W_{in}$, venant changer l'échelle des données en entrée.

In [55]:
reservoir = reset_esn()

states = []
scalings = [0.01, 0.1, 1.]
for iss in scalings:
    Win = mat_gen.generate_input_weights(units, 1, input_scaling=iss,
                                     proba=input_connectivity, input_bias=True,
                                     seed=seed)
    reservoir.Win = Win
    s = reservoir.compute_all_states([X_test1[:500]])
    states.append(s[0])
In [57]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(scalings)*100+10+i+1)
    plt.plot(s[:, :units_nb], alpha=0.6)
    plt.ylabel(f"$iss={scalings[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()

Correlation moyenne des activités des neurones du reservoir avec les entrées :

In [58]:
for i, s in enumerate(states):
    corr = correlation(states[i], X_test1[:500])
    print(f"ISS : {scalings[i]}, correlation moyenne : {corr}")

ISS : 0.01, correlation moyenne : -0.032300266438732315
ISS : 0.1, correlation moyenne : 0.5455501730434691
ISS : 1.0, correlation moyenne : 2.441693673875052
  • $+$ input scaling $\rightarrow$ activités corrélées aux données
  • $-$ input scaling $\rightarrow$ activités libres

L'input scaling peut aussi être utilisé pour ajuster l'influence de chaque donnée en entrée.

2.3. Le taux de fuite de la mémoire (leaking rate)¶

$$ x(t+1) = \underbrace{\color{red}{(1 - \alpha)} x(t)}_{\text{état actuel}} + \underbrace{\color{red}\alpha f(u(t+1), x(t))}_{\text{données suivantes}} $$

avec $\alpha \in [0, 1]$ et:

$$ f(u, x) = \tanh(W_{in} \cdotp u + W \cdotp x) $$
In [59]:
reservoir = reset_esn()

Win = mat_gen.generate_input_weights(units, 1, input_scaling=0.1, input_bias=True)
reservoir.Win = Win

states = []
rates = [0.02, 0.2, 0.9]
for lr in rates:
    reservoir.lr = lr
    s = reservoir.compute_all_states([X_test1[:500]])
    states.append(s[0])
In [60]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(rates)*100+10+i+1)
    plt.plot(s[:, :units_nb] + 2*i)
    plt.ylabel(f"$lr={rates[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()
  • $+$ leaking rate $\rightarrow$ faible inertie, faible mémorisation des états précédents
  • $-$ leaking rate $\rightarrow$ forte inertie, grande mémorisation des états précédents

Le leaking rate contrôle la "mémoire" de l'ESN. Il peut être vu comme l'inverse de sa constante de temps